HTIC  FILE  COPY  AD-A159  595 


► 

i 


WASHINGTON 
UNIVERSITY 
IN  ST  LOUIS 


REPORT  WU/CCM-84/1 


Contract  No.  N00014-81-K.0625 


Estimation  and 

Control  of  Error  Based  on 

P-Convergence 


B.  A.  Szabo 


JUNE,  1984 


CENTER  FOR 
COMPUTATIONAL 
MECHANICS 


DT®, 

^Et-ECTE 

W  sEPa'J's® 


WASHINGTON  UNIVERSITY 
CAMPUS  BOX  1129 
ST.  LOUIS,  MO  63130 


85  9  09  131 


Report  No.  WU/CCM-84/1 


Estimation  and  Control  of  Error 
Based  on  P-Convergence 


Barna  A.  Szabo 


June,  1984 


Accession  For 

KTIS 

GRA&I  M 

DTIC 

TAB  ^3 

Unannounced  Q 

i  c  a  t  i  1 

r 

By _ 

Distr 

Aval 

I’-'!  ■■  ■  *  y  'odes 

Avr, .  or 

Dlst 

id 

! 

1 

> 

1 

Paper  presented  at  the  International  Conference  on  Accuracy  Estimates  and 
Adaptive  Refinements  in  Finite  Element  Computations  (ARFEC)  Lisbon,  Portugal 
19-22  June,  1984. 


TABLE  OF^  CONTENTS  ; 


Summary 

1.  Introduction  . 

2.  Convergence  in  energy  norm' . 

3.  ^Problems  with  smooth  solutions' . 

'  3.1  Thick  walled  cylinder  under  internal  pressure 

3.2  Circular  hole  in  a  rectangular  panel  . 

4.  Problems  with  stress  singularities^ j . . . 

4.1  Short  cantilever  beam . .J. . 

4.2  Cracked  panel  . 

5.1-  Estimates  of  error  in  energy  norm  . 

5.1  Thick  walled  cylinder  under  Internal  pressure 

5.2  Cracked  panel  . 

6.  Summary  and  Conclusions  . 

7.  Acknowledgements  . . . . . 


8.  References 


ESTIMATION  AND  CONTROL  OF  ERROR  BASED  ON  P-CONVERGENCE 


SUMMARY 


The  relative  error  in  energy  can  be  estimated  on  the  basis  of  the 
asymptotic  estimate  of  the  rate  of  p-convergence.  Nearly  optimal  control 
of  error  can  be  achieved  when  feedback  information,  generated  from 
p-convergence  data,  is  used  for  deciding  whether  the  mesh  should  be 


refined  or  the  polynomial  degree  of  elements  increased.  -  , 


1.  INTRODUCTION 

In  order  to  keep  the  discussion  in  focus,  we  shall  be  concerned 
only  with  problems  of  plane  elasticity  and  finite  element  models  based 
on  the  displacement  formulation.  In  a  number  of  cases  the  domain  has 
reentrant  corners  or  sudden  changes  occur  in  the  boundary  conditions. 
In  the  neighborhood  of  such  points  the  exact  solution  is  of  the  form: 


u  =  I  A  r  ^  F  (0),  a  >  0 
i=l 


(1) 


where  u  is  the  displacement  vector,  r  and  9  are  polar  coordinates 
centered  on  the  point,  are  determined  from  the  condition  that  the 
solution  must  satisfy  the  Navier-Lame  equations  and  the  boundary  condi¬ 
tions  on  the  edges  that  meet  at  the  corner  [1].  For  the  purposes  of  the 
following  discussion  we  shall  assume  that  have  been  ordered  such  that 
a,  <  a„  <  a_  etc.  F.  are  sinusoidal  (vector)  functions.  A.  are  unknown 
coefficients  (amplitudes),  closely  related  to  the  stress  intensity 
factors  in  linear  elastic  fracture  mechanics.  In  fact,  in  linear  elastic 
fracture  mechanics  ~  “2  *  Y  ^1  *  »  ^2  * 

Kjj  being,  respectively,  the  Mode  I,  Mode  II  stress  intensity  factors. 

We  shall  be  interested  in  the  lowest  value  of  a^,  specifically  a^. 

Some  typical  values  that  frequently  occur  in  engineering  practice  are 
shown  in  Fig.  1.  When  is  large,  or  when  are  integers,  or 

when  the  points  of  singularity  lie  outside  of  and  far  from  the  solution 
domain  then  the  solution  is  said  to  be  smooth. 


Elastic  Solid 


a, «  0.6157 


Fig.  1 

Some  values  for  isotropic  elastic  solids. 

The  finite  element  solution  u  minimizes  the  potential  energy 

— ^FE 

expression  with  respect  to  a  set  of  functions  that  can  be  written  in  the 
following  form: 

N 

u  ”  2!  (2) 

i=l 

where  are  the  basis  functions,  constructed  from  the  element 

shape  functions  in  such  a  way  that  the  appropriate  continuity  and  boundary 


conditions  are  satisfied.  The  basis  functions  depend  on  the  choice  of 
the  finite  element  mesh  for  the  solution  domain,  the  shape  functions 
defined  on  the  finite  elements  and  the  mapping  of  the  elements: 

X  “  x(5,  n),  y  =  y(5.  n)  (3) 

where  5  b  are  the  standard  coordinates  and  (3)  represents  transforma¬ 
tion  of  the  standard  element  into  the  'real'  element.  The  element  shape 
functions  are  defined  on  the  standard  element.  We  shall  be  concerned 
only  with  polynomial  shape  functions.  Unless  the  transformation  is 
isoparametric  or  subparametric  however,  the  mapped  shape  functions  and 
therefore  the  basis  functions  are  not  polynomials.  The  set  of  all 
functions  that  can  be  written  in  the  form  of  eq.  (2)  is  called  a  finite 
element  space. 

A  particular  choice  of  mesh  and  polynomial  degree  p  may  not  yield 
an  approximate  solution  of  sufficient  precision.  It  is  necessary  there¬ 
fore  to  have  the  ability  to  improve  the  quality  of  approximation.  This 
involves  increasing  the  number  of  degrees  of  freedom.  When  the  number 
of  degrees  of  freedom  is  increased  in  such  a  way  that  the  finite  element 
spaces  of  fewer  degrees  of  freedom  are  embedded  in  the  finite  element 
spaces  of  greater  number  of  degrees  of  freedom  then  the  process  of 
increasing  the  degrees  of  freedom  is  called  extension.  The  algorithm  or 
strategy  by  which  the  number  of  degrees  of  freedom  is  increased  is 
called  extension  operator. 

Extension  operators  are  classified  into  three  major  categories 
called  h-extension,  p-extension,  and  h-p  extension.  When  aspects  of 
implementation  are  emphasized  then  the  word  version  rather  than  extension 
is  used.  In  the  h-version  the  polynomial  degree  of  the  elements  is 
fixed,  therefore  extension  is  possible  only  by  mesh  refinement.  The 
size  of  the  elements  is  usually  denoted  by  h,  hence  the  name:  h-version. 
In  the  p-version  the  polynomial  degree  of  the  elements  may  be  varied 
over  a  substantial  range.  Therefore  extension  is  possible  through 
Increasing  the  polynomial  degree  of  elements.  In  the  h-p  extension  the 
number  of  degrees  of  freedom  is  increased  by  some  combination  of  mesh 
refinement  and  increase  of  polynomial  degree. 


-4- 


2.  CONVERGENCE  IN  ENERGY  NORM 

Finite  element  solutions  minimize  the  error  in  energy  norm.  There¬ 
fore  the  magnitude  of  the  error  in  energy  norm  is  a  good  measure  of  the 
overall  quality  of  approximation.  Asymptotic  error  estimates  are  avail¬ 
able  for  the  three  methods  of  extension.  Following  is  a  brief  summary: 

In  the  h-extenslon  when  a  sequence  of  meshes  is  obtained  through 
uniform  refinement  then  the  estimate  of  error  is: 

^1 

-  jjl/2  min(p,a) 

where  ] |e| [j,  is  the  error  in  energy  norm  (i.e.  the  square  root  of  the 
strain  energy  of  the  error);  N  is  the  number  of  degrees  of  freedom;  p  is 
the  polynomial  degree  of  elements  and  a  measures  the  smoothness  of  the 
exact  solution.  Specifically,  in  the  problems  considered  herein,  a  is 
the  smallest  exponent  of  r  in  eq.  (1),  i.e.  a  =  a^.  The  exponent  of  N 
is  called  the  asymptotic  rate  of  convergence. 

When  the  sequence  of  meshes  is  obtained  in  such  a  way  that  the 
error  associated  with  each  element  is  approximately  the  same  then  the 
meshes  are  called  equilibrated  meshes  and  the  estimate  is: 


(5) 


Note  that  the  rate  of  convergence  is  independent  of  ct  and  can  be  made 
arbitrarily  large  by  choosing  sufficiently  large  p. 

In  p-extension  the  mesh  is  fixed  and  the  polynomial  degree  of 
elements  is  Increased.  In  this  case,  when  the  point  of  singularity  is 
located  at  a  nodal  point,  the  estimate  is: 


E 


(6) 
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Note  that  the  rate  of  convergence  is  exactly  twice  the  rate  of  convergence 
of  the  h-extension  based  on  uniform  mesh  refinement. 

In  the  h-p  extension  optimal  mesh  refinement  is  coupled  with  optimal 
p-distribution  and  the  estimate  is 


0 

expCyN  ) 


(7) 


where  y ,  0  are  constants,  0  >  1/3.  In  this  case  the  asymptotic  rate  of 
convergence  is  exponential. 

Some  of  these  results  are  illustrated  graphically  in  Fig.  2.  The 
relative  errors  shown  in  Fig.  2  are  typical  values. 


C 

a> 

u 

w 

(U 

Q. 

(T 

O 

ac 

cc 

Ui 

UJ 

> 

UJ 

(T 


Fig.  2 

Performance  of  the  h-,  p-  and  h-p  extension  processes. 

(The  relative  error  values  shown  are  typical  for  certain  engineering  problems) . 
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The  asymptotic  rate  of  convergence  of  p-extensions  was  discussed  in 
a  number  of  papers  and  several  examples  were  presented  [2-7].  In  [3]  a 
rigorous  proof  of  eq.  (6)  was  given.  In  the  examples  the  meshes  were 
usually  designed  in  such  a  way  that  the  minimum  number  of  elements 
needed  for  representing  the  domain  were  used.  An  important  exception 
was  a  demonstration  of  the  h-p  extension  in  [4]. 

In  this  paper  it  will  be  shown  that:  (1)  the  exponential  rate  of 
convergence  is  in  fact  realized  in  the  case  of  problems  with  smooth 
solutions  when  the  p-extension  is  used;  (2)  Problems  with  stress  singu¬ 
larities  can  be  made  to  behave  almost  as  problems  with  smooth  solutions 
do  when  properly  graded  meshes  are  used  in  conjunction  with  the  p-extension 
and  (3)  the  relative  error  in  energy  norm  can  be  estimated  on  the  basis 
of  the  asymptotic  estimate  for  the  p-extension,  eq.  (6). 


3. 


PROBLEMS  WITH  SMOOTH  SOLUTIONS 


3.1  Thick  walled  cylinder  under  internal  pressure. 

Let  us  first  consider  a  thick  walled  cylinder  under  internal  pressure. 

In  this  case  the  classical  solution  is  known.  (See,  for  example,  [8]). 

Assuming  plane  strain  conditions  and  Poisson's  ratio  of  0.3,  the  strain 

2  2 

energy  U  for  a  45  degree  sector,  as  shown  in  Fig.  3,  is  0.748746249  p  r^  /E 
per  unit  length  of  the  cylinder,  given  of  course  that  r^/r^  =  2. 

E  is  the  modulus  of  elasticity. 

The  relative  error  in  energy  norm  is: 


(8) 


where:  U  exact  strain  energy 

U  computed  strain  energy,  polynomial  degree  p. 


Fig.  3 

45  degree  sector  of  a  thick  walled 
cylinder  under  internal  pressure. 
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Using  a  single  quadrilateral  element,  mapped  by  linear  blending  so  that 
the  annular  sector  is  represented  exactly,  (e  )_  was  computed  for 
p  =  1,2,... 8.  The  results  are  shown  in  Figs.  4  and  5.  In  Fig.  4  the 
relative  error  is  plotted  against  the  number  of  degrees  of  freedrn,  N, 
on  log-log  scale.  It  is  seen  that  the  slope  of  the  curve  increases  with 
N.  This  is  consistent  with  the  qualitative  illustration  of  convergence 
in  Fig.  2.  In  Fig.  5  the  relative  error  is  plotted  against  the  number 
of  degrees  of  freedom  on  semi-log  scale.  In  this  case  the  slope  of  the 
relative  error  approaches  a  constant  value,  indicating  that  the  rate  of 
convergence  is  nearly  exponential. 

Note  that  p  =  4  is  more  than  sufficient  for  engineering  accuracy 
for  this  problem. 


NUMBER  OF  DEGREES  OF  FREEDOM 


Fig.  4 

45  degree  sector  of  a  thick  walled  cylinder  under  internal  pressure. 
Relative  error  in  energy  norm  plotted  against  N  on  log-log  scale. 


NUMBER  OF  DEGREES  OF  FREEDOM 


Fig.  5 

45  degree  sector  of  a  thick  walled  cylinder  under  internal  pressure. 
Relative  error  in  energy  norm  plotted  against  N  on  semi-log  scale. 


3.2  Circular  hole  in  a  rectangular  panel 

Let  us  next  consider  the  well  known  problem  of  a  circular  hole  in  a 
rectangular  panel  subjected  to  uniaxial  tension.  See,  for  example  [9]. 
Because  the  solution  of  this  problem  is  smooth  (no  reentrant  corners  are 
present  and  no  sudden  changes  in  material  properties  or  loading  conditions 
occur)  the  minimum  number  of  elements  needed  for  representing  the  geometry 
is  the  best  choice  for  the  p-version.  In  the  h-version  the  error  is  con¬ 
trolled  by  mesh  refinement.  A  typical  h-version  mesh  (produced  by  a 
practicing  engineer)  and  the  p-version  mesh  are  shown  in  Fig.  6.  Of 
interest  is  the  maximum  stress.  The  maximum  stress,  computed  directly 
from  the  finite  element  solution  using  the  three  element  mesh  with  the 
polynomial  degree  ranging  from  1  to  8  and  —  =  0.5,  is  shown  in  Fig.  7. 
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6.  SUMMARY  AND  CONCLUSIONS 

When  properly  graded  meshes  are  used  then  the  performance  of  the 
p-extension  in  the  pre-asymptotic  range  is  very  close  to  the  best  per¬ 
formance  attainable  by  sequences  of  optimal  meshes  coupled  with  sequences 
of  optimal  p-distributions. 

The  p-extension  process  provides  means  for  estimating  the  relative 
error  in  energy  norm.  It  also  provides  information  on  the  basis  of 
which  it  is  possible  to  decide  whether  the  mesh  should  be  refined  or  the 
polynomial  degree  of  elements  should  be  increased. 

When  the  exact  solution  is  smooth  then  the  coarsest  possible  mesh 
should  be  used.  In  this  case  mesh  design  is  controlled  entirely  by  the 
geometry  of  the  domain.  The  only  restriction  on  the  mesh  is  that  the 
mapping  of  the  elements  must  be  smooth  also,  therefore  large  aspect 
ratios  must  be  avoided.  In  the  p-extensions  aspect  ratios  as  large  as  20:1 
are  permissible,  however. 

When  the  exact  solution  is  not  smooth,  for  example  corner  singular¬ 
ities  are  present,  then  the  points  of  singularity  must  be  isolated  by 
one  or  more  layers  of  small  elements.  Grading  of  the  elements  should  be 
such  that  the  element  sizes  are  in  geometric  progression  with  the  smallest 
element (s)  at  the  point  of  singularity.  The  geometric  progression 
should  have  a  common  factor  of  about  0.15.  In  this  case  entry  into  the 
asymptotic  range  is  shifted  toward  higher  p-values  and  the  pre-asymptotic 
convergence  of  the  p-extension  process  is  much  stronger  than  asymptotic. 

When  the  mesh  is  properly  designed  then  the  finite  element  solution 
behaves  almost  as  if  the  exact  solution  were  smooth. 

In  general,  overref inement  in  the  neighborhood  of  singular  points 
is  preferable  to  underref inement. 

Proper  mesh  design  is  highly  beneficial  from  the  point  of  view  of 
stress  computations:  The  oscillatory  behavior  of  stresses  can  be  attenuated 
to  such  an  extent  that,  with  the  exception  of  the  immediate  neighborhood  of 
singular  points,  the  error  is  very  small  everywhere. 
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TABLE  IV 

Estimated  and  true  error  in  energy  norm. 
Cracked  panel.  Mesh  B. 


£ 

N 

_£ 

U  E 

e 

Percent  Rel. 
Est.  'd 

Error 

True 

1 

27 

0.207859 

_ 

. 

35.10 

2 

78 

0.2330A8 

- 

- 

13.02 

3 

137 

0.236076 

1.71 

8.89 

6.46 

A 

220 

0.236739 

2.55 

3.A6 

3.71 

5 

327 

0.236895 

2.87 

1.76 

2.68 

6 

A58 

0.236953 

2.22 

1.48 

2.17 

7 

613 

0.236983 

1.66 

1.42 

1.86 

8 

792 

0.237001 

1.26 

1.43 

1.64 

Note  that  the  preasymptotic  range  is  extended  when  Mesh  B  is  used. 

Systematic  application  of  the  procedure  just  described  is  the 
'feedback  h-p  extension'  in  the  terminology  of  reference  [12]. 
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It  is  seen  that  g  is  monotonically . increasing.  Although  the 
relative  errors  in  energy  norm  are  very  small,  the  estimated  and  true 
errors  are  close. 


5.2  Cracked  panel 

We  first  tabulate  the  convergence  data  for  Mesh  A  (see  Table  III). 
The  observed  rate  of  convergence  (6)  is  clearly  stronger  than  the 
asymptotic  rate  (2a^  =  1)  up  to  p  =  6.  For  higher  p  values  the  rate 
of  convergence  becomes  weaker  and  approaches  the  asymptotic  rate.  Note 
also  that  the  estimated  and  true  error  values  are  reasonably  close. 

At  p  =  8  the  relative  error  in  energy  norm  is  4.22  percent  relative 
error  in  strain  energy.  From  the  practical  point  of  view  this  level  of 
precision  is  usually  more  than  sufficient.  For  example,  the  relative 
error  in  the  stress  intensity  factor  K^,  computed  from  the  strain  energy 
release  rate  at  p  =  8,  is  only  0.38  percent.  Nevertheless,  if  greater 


TABLE  III 


Estimated  and  true  error  in  energy  norm. 
Cracked  panel.  Mesh  A. 


£ 

N 

_E 

U  E 
_E_ 

1 

18 

0.201706 

2 

52 

0.229199 

3 

94 

0.234117 

4 

152 

0.235489 

5 

226 

0.236053 

6 

316 

0.236346 

7 

422 

0.236525 

8 

544 

0.236643 

Percent  Rel.  Error 


3 

Est.  'd 

True 

38.62 

- 

- 

18.21 

1.31 

13.38 

11.15 

1.96 

6.11 

8.15 

1.57 

5.26 

6.53 

1.33 

4.70 

5.51 

1.10 

4.49 

4.77 

1.02 

4.11 

4.22 

precision  is  needed,  then  the  mesh  should  be  refined  using  the  geometric 
progression  rule.  Repeating  the  procedure  for  Mesh  B  shown  in  Fig.  11, 
the  results  shown  in  Table  IV  were  obtained. 
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the  h-extension  process.  In  fact,  the  relative  error  for  the  cantilever 
problem  shown  in  Fig.  10  was  estimated  in  this  way.  When  p-extension  is 
used  then  B  When  h-extension  is  used  with  uniform  mesh  refinements 

then  6  =  min(a^,p). 

Referring  to  Fig.  2,  we  note  that  when  the  mesh  is  well  designed 
then  the  extension  is  not  in  the  asymptotic  range.  Therefore  the  meshes 
should  be  designed  in  such  a  way  that  the  asymptotic  range  is  not  entered 
when  the  polynomial  degree  is  increased.  Ideally,  the  mesh  should  be 
designed  so  as  to  cause  6  to  Increase  with  N.  When  three  successive 
finite  element  solutions  are  available  in  the  extension  process,  i.e. 

U  ,  U  I ,  U  and  N  ,  N  , ,  N  „  are  known,  then  we  can  estimate  3  by 

p  p-1  p-2  p  p-1  p-2 

assuming  C  to  be  constant  and  solving  the  following  equation  for  3: 


U  - 
_P _ 


U  1  -  - 

P-J-  P-2 

N-^  -  N-^ 

p-2  p-1 


(14) 


If  we  find  that  3  is  not  increasing  with  p  and  the  error  is  still  large 
then  those  elements  that  have  vertices  on  the  singular  point  should  be 
refined,  using  the  geometric  progression  rule  for  refinement.  The 
procedure  is  illustrated  on  the  basis  of  two  problems:  one  has  a  very 
smooth  solution,  the  other  is  characterized  by  a  strong  stress  singularity. 

5.1  Thick  walled  cylinder  under  internal  pressure. 

The  problem  was  solved  with  a  single  finite  element  which  is  shown 
in  Fig.  3.  The  convergence  data  are  listed  in  Table  II. 


TABLE  II 


Estimated  and  true  error  in  energy  norm. 
Thick  walled  cylinder  under  internal  pressure. 


U  E 


Percent  Rel.  Error 


£ 

N 

_E 

3 

Est. 'd 

True 

1 

4 

0.6539035818 

35.5905 

2 

10 

0.7421867856 

- 

- 

9.3598 

3 

16 

0.7483608276 

2.6312 

5.8098 

2.2688 

4 

24 

0.7487274398 

5.9418 

0.6954 

0.5012 

5 

34 

0.7487454062 

7.3680 

0.1413 

0.1061 

6 

46 

0.7487462140 

8.8350 

0.0283 

0.0217 

7 

60 

0.7487462477 

10.4389 

0.0055 

0.0043 

8 

76 

0.7487462491 

12.0693 

0.0010 

0.0008 
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5.  ESTIMATES  OF  ERROR  IN  ENERGY  NORM 

The  p-extension  process  can  be  used  for  obtaining  valuable  (feedback) 
Information  concerning  the  overall  quality  of  approximation.  On  the 
basis  of  this  Information  It  Is  possible  to  estimate  the  relative  error 
In  energy  norm  and  decide  whether  It  Is  better  to  refine  the  mesh  or 
Increase  the  polynomial  degree  of  elements. 

The  asymptotic  estimates  of  error  In  energy  norm  were  reviewed  In 
Section  2.  Equivalently,  the  error  In  strain  energy  can  be  written  as: 


|U-U  1  <  cn"® 
I  pi  _  p 


(10) 


where  U  Is  the  exact  strain  energy,  Is  the  computed  strain  energy,  N^ 

Is  the  corresponding  number  of  degrees  of  freedom,  subscript  p  represents 

the  pol5m.omlal  degree  of  elements,  C  Is  a  positive  constant  and  3  Is  the 

rate  of  convergence.  In  the  case  of  homogeneous  kinematic  boundary 

conditions  U  >  U  and  we  can  write: 

P 


u-u  <  cn"^. 

p  -  p 


(11) 


When  the  finite  element  solution  Is  In  the  asymptotic  range  and  3  Is 
known  then  we  can  readily  compute  C  from: 


C 


and  therefore  U  can  be  estimated  from: 


U 


UN^ 
P  ,P 


-  U  ,N 

P-1 


-  N 


P-1 


8 

P^l 


(12) 


(13) 


Although  the  exact  solution  Is  not  known,  the  strain  energy  of  the  exact 
solution  can  be  estimated  with  great  precision  through  either  the  p-  or 
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Since  the  exact  solution  is  known,  we  are  in  position  to  examine  the 
accuracy  of  stress  computations.  Ue  have  computed  the  stresses  at 
X  =  0.25a  y  *  0.25a  (See  Fig.  11).  Because  this  point  lies  at  the 
boundary  of  two  elements,  we  obtain  two  sets  of  stress  values,  depending 
on  which  element  is  used  for  the  stress  computations.  We  refer  the 
element  with  vertices  at  (0,0),  (a, a)  (0,a)  as  element  1,  the  element 
with  vertices  at  (0,0)  (a,0),  (a, a)  as  element  2.  The  computed 
values  for  both  elements,  together  with  the  corresponding  relative 
errors,  are  shown  in  Table  I.  Note  that  the  error  is  very  small  for  a 
wide  range  of  p-values. 

TABLE  I 

Stress  Oy  computed  at  x  =  0.25a,  y  =  0.25a, 


using 

Mesh  A.  The  exact  value 

of  o  is: 

y 

0.8390  Kj/Ja 

Element  1 

Element  2 

£ 

0  \[a/K- 
y^  I 

Rel.  Err.  (%) 

a 

y^  I 

Rel .  Err . 

1 

0.7891 

5.95 

0. 5804 

30.8 

2 

0.7336 

12.6 

0.7516 

10.4 

3 

0.8355 

0.42 

0.8510 

1.43 

4 

0.8456 

0.79 

0.8411 

0.25 

5 

0.8402 

0.14 

0.8407 

0.20 

6 

0.8359 

0.37 

0.8383 

0.08 

7 

0.8352 

0.45 

0.8362 

0.33 

8 

0.8359 

0.37 

0.8361 

0.35 

Alternative  means  of  stress  computation  are  available  [10,11].  The 
results  shown  here  indicate  however  that  direct  computation  of  stresses 
from  finite  element  solutions  is  adequate  for  practical  purposes  when 
the  mesh  is  properly  designed.  This  is  especially  important  in  the  case 
of  design  computations  where  the  overall  behavior  of  the  solution  is  of 
Interest.  The  extraction  techniques  presented  in  [10,11]  and  investigated 
in  [6,7]  are  of  great  practical  Importance  when  the  solution  in  the 
immediate  neighborhood  of  a  stress  singularity  is  of  interest. 


RELATIVE  ERROR  IN  ENERGY  NORM 


•nr 


cos  sin  sin  -j J  ,  -n  _<  6  _<  n 


(9b) 


K, 


xy 


nr 


6  0  36 

sin  cos  cos  ~2  *  -tt  £  6  £  n. 


(9c) 


Therefore  the  exact  solution  is  known  and  the  strain  energy  cam  be 

2 

computed.  It  is  0.23706469  Kj.a/E  for  the  half  panel.  (Due  to  symmetry 
only  half  of  the  panel  has  to  be  ainalyzed). 

The  p-convergence  paths  for  Meshes  A  and  B  are  shown  in  Fig.  12. 
Note  that  the  convergence  path  corresponding  to  Mesh  B  is  very  similar 
to  the  qualitative  behavior  shown  in  Fig.  2  for  strongly  graded  meshes. 


10  100  1000 
NUMBER  OF  DEGREES  OF  FREEDOM 

Fig.  12 

Cracked  panel.  Performance  of  the  p-extension 
when  strongly  graded  meshes  are  used. 


in  the  number  of  degrees  of  freedom  when  two  elements  are  added  is  not 
substantial  (Fig.  10).  Similarly,  the  use  of  graded  p-distribution, 
i.e.  assigning  low  polynomial  degrees  to  the  small  elements  at  the 
singularities  and  high  polynomial  degrees  to  the  large  elements  away 
from  the  singularities,  as  in  the  h-p  extension,  would  not  produce 
appreciable  savings  in  practical  computations. 

4 . 2  Cracked  Panel 

Let  us  now  consider  the  cracked  panel  shown  in  Fig.  11. 


0.0225a 


MESH  A 


MESHB 


Fig.  11 
Cracked  panel. 

We  assume  plane  strain  conditions,  Poisson's  ratio  of  0.3,  and  unit 
thickness  and  load  the  panel  in  such  a  way  that  the  tractions  exactly 
correspond  to  the  well  known  first  term  of  the  asymptotic  expansion  for 
Mode  1 ,  i.e.: 


9  n  4  e  J  30, 

cos  TT  [1  -  sin  rr  sin  — ;r]  .  -ir  <  6  <  ir 
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The  rates  of  convergence  in  energy  norm  are  shown  in  Fig.  10  for 
the  h-extension  with  uniform  mesh  refinement  (p=l) ,  and  the  p-extension 
for  three  different  mesh  layouts.  It  is  seen  that  the  theoretically 
predicted  rate  is  realized  in  the  case  of  both  the  h-  and  p-versions 
when  the  mesh  is  uniform,  i.e.  Mesh  A  is  used.  In  either  case  it  is 
impractical  to  achieve  1  percent  relative  error  in  energy  norm.  The 
rate  of  convergence  is  so  small  that  the  h-  and  p-extension  processes 
do  not  provide  effective  control  of  error  in  this  norm.  On  the  other 
hand,  the  convergence  of  the  p-extension  process  is  very  strong  when 
strongly  graded  meshes  are  used.  Mesh  B5  is  a  five-element  mesh,  shown 
in  Fig.  9,  Mesh  B7  is  a  seven-element  mesh,  the  grading  of  the  elements 
is  in  geometric  progression  toward  the  singularities,  the  common  factor 
being  0.1.  The  strong  convergence  obtained  with  these  meshes  seemingly 
contradicts  the  theoretically  predicted  rate  of  eq.  (6).  This  can  be 
explained  as  follows:  The  estimates  given  by  eqs.  (4)  to  (7)  are 
asymptotic  estimates,  hence  they  are  valid  only  for  'large'  N  values. 

The  size  of  N  for  the  estimates  to  hold  depends  on  the  mesh  design. 

When  the  elements  at  singular  points  are  large  then  the  error  of  approxi¬ 
mation  is  controlled  by  the  error  associated  with  these  elements  and  N 
does  not  have  to  be  very  large  for  the  asymptotic  estimates  to  govern. 
When  the  elements  at  the  singular  points  are  small,  as  in  the  case  of 
Mesh  B5,  the  error  associated  with  the  larger  elements  away  from  the 
singularity,  where  the  solution  is  smooth,  is  the  controlling  factor. 
Consequently  the  finite  element  solution  behaves  as  if  the  exact 
solution  were  smooth.  Entry  into  the  asymptotic  range  occurs  only 
when  the  polynomial  degree  is  sufficiently  high  so  that  the  elements 
at  singular  points  begin  to  control  the  error  of  approximation.  The 
validity  of  this  explanation  is  confirmed  by  the  results  obtained  with 
Mesh  B7:  Refinement  at  the  singularity  does  not  significantly  reduce 
the  error.  In  fact,  for  low  p-values  it  is  more  efficient  to  use 
Mesh  B5  than  Mesh  B7.  The  cross-over  occurs  at  about  p  “  5. 

The  strong  preasymptotic  behavior  of  the  p-extension  associated 
with  Meshes  B5,  B7  is  similar  to  the  asymptotic  behavior  of  problems 
with  smooth  solutions. 

Note  that  overrefinement  (Mesh  B7)  is  not  detrimental  from  the 
point  of  view  of  roundoff  error  and  is  not  very  wasteful.  The  Increase 


LOG(er) 


u 

L 


-13- 
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Fig.  9 

Short  cantilever  beam. 
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4.  PROBLEMS  WITH  STRESS  SINGULARITIES 

When  the  mlnlmuin  number  of  elements  Is  used  and  stress  singularities 
are  present  (i.e.  the  form  of  the  exact  solution  at  one  or  more  points 
is  that  of  eq.  (1)  with  less  than  one)  then  entry  into  the  asymptotic 
range  occurs  at  low  p-values,  usually  p  =  3  or  p  =  4,  and  the  stress 
field  exhibits  strong  oscillatory  behavior.  This  oscillatory  behavior 
is  closely  related  to  the  smoothness  of  the  exact  solution:  The  smoother 
the  exact  solution,  the  less  pronounced  are  the  oscillations.  These 
oscillations  are  beneficial  in  the  sense  that  the  faster  rate  of  con¬ 
vergence  of  the  p-extension  in  energy  (as  compared  with  the  h-extension 
based  on  uniform  meshes)  is  owed  to  the  property  of  polynomials  that 
they  are  able  to  oscillate  with  increased  frequency  as  the  polynomial 
degree  is  increased  and  the  wavelength  of  oscillations  decreases  with 
distance  measured  from  element  boundaries.  On  the  other  hand,  stress 
oscillations  are  confusing  when  stresses  are  of  primary  interest. 

Numerical  experiments  have  shown,  however,  that  the  boundaries  of  finite 
elements  attenuate  stress  oscillations.  The  attenuation  is  so  substantial 
that  with  proper  mesh  design  the  error  in  stress  maxima,  outside  of  the 
immediate  neighborhood  of  stress  singularities,  can  be  reduced  to  under 
1  percent  for  all  stress  analysis  problems  of  practical  Importance. 

In  the  h-p  extension  the  meshes  are  strongly  graded  toward  the 
singularity;  the  element  sizes  are  reduced  in  geometric  progression  with 
a  common  factor  of  about  0.15  and  the  polynomial  degrees  of  elements  are 
assigned  such  that  the  smallest  element  has  the  lowest  polynomial 
degree,  the  largest  element  the  highest.  Nearly  as  good  performance  can 
be  achieved  when  the  elements  are  strongly  graded  toward  the  singularity 
but  the  distribution  of  polynomial  degrees  is  uniform.  In  the  following 
we  demonstrate  the  performance  of  the  p-extension  on  the  basis  of  two 
examples . 

4.1  Short  Cantilever  Beam 

The  geometric  definition  of  the  problem  and  the  meshes  used  ir  our 
analysis  are  shown  in  Fig.  9.  Plane  strain  conditions  and  Poisson's 
ratio  of  0.3  are  assumed. 
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lt  is  seen  that  there  are  minor  oscillations  in  the  computed  stresses 
but  for  p  *  4  to  8  the  stress  value  is  well  within  the  5  percent  relative 
error  range.  For  p  =  8  the  finite  element  solution  is  identical  with 
the  solution  presented  in  [9].  More  importantly,  the  extension  process 
confirms  the  fact  that  convergence  has  occurred.  The  variation  in 
stresses  is  0.3  percent  in  the  range  of  p  =  6,7,8.  Given  the  fact  that 
the  solution  is  smooth  and  therefore  convergence  (in  energy)  is  strong, 
one  can  reasonably  estimate  that  the  relative  error  is  under  1  percent. 

The  fact  that  the  p-version  tolerates  large  aspect  ratios  is  illus¬ 
trated  in  Fig.  8.  Using  the  same  three  element  mesh  as  shown  in  Fig.  6, 
the  stress  concentration  factor  was  computed  for  a  wide  range  of  r/w 
ratios.  It  is  seen  that  substantial  deviation  from  the  theoretical 
value  occurs  only  at  very  large  and  very  small  r/w  ratios.  The  aspect 
ratios  of  the  finite  elements  are  very  large  in  those  cases,  however. 


Fig.  8 

r 

Stress  concentration  factor  as  function  of  — . 
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Fig.  6 

Quarter  model  of  a  rectangular  panel  with  a  circular  hole 


Fig.  7 

Ratio  of  maximum  stress  to  applied  stress 
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